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We present an implementation of the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional within the 
full-potential linearized augmented-plane-wave (FLAPW) method. Pivotal to the HSE functional is 
the screened electron-electron interaction, which we separate into the bare Coulomb interaction and 
the remainder, a slowly varying function in real space. Both terms give rise to exchange potentials, 
which sum up to the screened nonlocal exchange potential of HSE. We evaluate the former with 
the help of an auxiliary basis, defined in such a way that the bare Coulomb matrix becomes sparse. 
The latter is computed in reciprocal space, exploiting its fast convergence behavior in reciprocal 
space. This approach is general and can be applied to a whole class of screened hybrid functionals. 
We obtain excellent agreement of band gaps and lattice constants for prototypical semiconductors 
and insulators with electronic-structure calculations using plane-wave or Gaussian basis sets. We 
apply the HSE hybrid functional to examine the ground-state properties of rocksalt GdN, which 
have been controversially discussed in literature. Our results indicate that there is a half-metal 
to insulator transition occurring between the theoretically optimized lattice constant at K and 
the experimental lattice constant at room temperature. Overall, we attain good agreement with 
experimental data for band transitions, magnetic moments, and the Curie temperature. 
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I. INTRODUCTION 

Density functional theory (DFTjPIl jg ^ powerful tool 
for calculating the electronic ground-state properties of 
molecules and solids. The predictive power of numerical 
DFT calculations relies on the availability of accurate 
approximations for the exchange-correlation (xc) energy 
Bxc, which incorporates all complicated many-body ef- 
fects. In many systems, this quantity is described ad- 
equately by the local-density approximation (LDA).'2HS1 
However, in more complex systems, physical properties 
such as the geometric structure, magnetic properties, and 
the band gap are not well reproduced. One can go beyond 
the LDA by taking into account the local density gradi- 
ent, which yields the generalized gradient approximation 
(GGA) ^ upon which many functionals are based. How- 
ever, despite their success, the GGA functionals still of- 
ten fail in describing systems with localized states, which 
is attributed to an incomplete cancellation of the self- 
interaction error in these semi-local functionals.'^ 

This deficiency is particularly critical in systems whose 
electronic properties are largely governed by the corre- 
lated motion of electrons in localized states. The rare- 
earth chalcogenides are among this class of materials, 
having iircompletely filled /-electron shells. They are iir- 
sulating, semiconducting, or metallic depending on de- 
tails of the valency of the rare-earth element. Gadolin- 
ium nitride (GdN) is widely studied owing to the ferro- 
magnetic order, large magnetic moment of 6.88ub per 
Gd atorcPSland its large magnetoresistive effect, El which 
makes the material interesting for technological applica- 
tions. The mechanism of the ferromagiretic order is still 
under debate. Variou s types are being discussed, such as 
carrier mediatecP^EH and superexchange mechanisms .1^ 
Another point of debate are the electronic properties. 



It was experimentally demonstrated to be a low car- 
rier semimetal in single crystali^ and insulating in thin 
films.'i^ There are also several recent reports of thin films 
of GdN having a degenerately doped semiconductiir^i^IIIl 
or a metallic ground stat^^^l based on the resistivity data 
measured at low temperatures. Theo retically it is pre- 
dicted t o have a semiconductin^^oiSU qj. half-metallic 
charactei'i^'^21111 based on ah initio calculations. 

Materials with strongly localized states, such as the 
/ states in GdN, are often treated within the LSDA-|-[/ 
methodjI^S where the electron correlation in these states 
is described with an additional on-site term, that in- 
volves the Hubbard parameter U . The disadvantage of 
this method is that the value of U is not known a pri- 
ori. Although methods to estimate the U value from 
first-priirciples calculations have been developed,'25H22l ]^ 
is usually chosen to reproduce experimental observations. 
However, a specific U value that provides a good descrip- 
tion of one quantity is often not suitable to describe an- 
other quantity.'^ 

During the last decade, hybrid functionals that com- 
bine a fraction of nonlocal Hartree-Fock (HF) exchange 
with local xc functionals have been shown to be a 
viable improvement over LDA and GGA offering a 
parameter-fre e des cription specifically suited for band 
gap materials.'^Slflni xhe explicit consideration of nonlo- 
cal HF exchange leads to a partial cancellation of the 
self-interaction error, but also makes numerical calcu- 
lations considerably more demanding than conventional 
LDA and GGA calculations. Various hybrid functionals 
have been developed. In empirical hybrid functionals, 
such as the B3LYP functional,'^^^ the fraction of HF ex- 
change is determined by fitting to an experimental data 
set. In the PBEO functionaP^ the mixing parameter for 
the HF exchange is inferred from expanding the inte- 
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grand of the adiabatic-connection formula of the exact 
xc functionaL 

In periodic systems, the Coulomb interaction between 
the electrons is effectively screened by polarization ef- 
fects in the electron system. The effective interaction is 
particularly short-ranged in systems with small or van- 
ishing band gaps. Therefore, Heyd et al^ introduced a 
range-separated hybrid functional, which has the added 
benefit to reduce the computational cost within a basis 
of localized Gaussian functions. Starting from the PBEO 
functional, they partitioned the Fock exchange term into 
a short- and a long-range part, where the former is de- 
scribed by a correspondingly screened Fock term and the 
latter is treated by a local approximation, derived from 
the PBE functional.'^Heyd et al. showed that this hybrid 
functional leads to a reduced computational demand for 
localized basis sets compared with the PBEO functional. 
Furthermore, it even yields results which are often in bet- 
ter agreement with experiment. 

The Heyd-Scuseria-Ernzerhof (HSE) hybrid functional 
has been implemented within Gaussian'^'^ and plane- 
wav^ basis sets. In this work, we present an imple- 
mentation within the full-potential linearized augmented- 
plane-wave (FLAPW) approach as implemented in the 
Fleur codCjISSl which provides a highly accurate all- 
electron basiJ^^MsH for a large variety of materials, includ- 
ing open systems with low symmetry, d- and /-electron 
systems, as well as oxides and nitrides. An implementa- 
tion of the PBEO functional limited to certain localized 
states and on-site interactions was given by Tran et al.^^ 
Betzinger et al.^^ described an efficient way to calculate 
the full nonlocal exchange potential for the PBEO func- 
tional without these restrictions. Very recently, Tran and 
Blah£p^ reported an implementation of hybrid function- 
als whose nonlocal exchange integrals are evaluated us- 
ing the pseudocharge method of Weinert.'^ However, this 
approach mathematically restricts the electron-electron 
interaction to potentials that are solutions of Laplace- 
type equations, whose radial solutions can be expressed 
as analytically or numerically known regular and irreg- 
ular solutions and spherical harmonics, such as the bare 
Coulomb and the screened Yukawa potential. The error 
function used in the HSE functional does not have this 
property. In our implementation, there is no such restric- 
tion. Our approach is very general. In fact, any kind of 
interaction potential can be implemented for the nonlo- 
cal exchange potential by changing only a single line of 
code. The only requirement is that it differs from the 
bare Coulomb potential by a function that possesses a 
fast Fourier expansion, a condition that is fulfilled by all 
physical screened potentials, including the error function 
used in the HSE functional. 

Our numerical approach extends the implementation of 
Ref. mi which is based on an auxiliary basis that is de- 
signed to represent products of wave functions. This so- 
called mixed product basis (MPB) is constructed directly 
from products of LAPW basis functions and retains the 
full accuracy of the all-electron description. Several tech- 



niques were introduced to accelerate the evaluation of the 
computationally expensive nonlocal exchange term. Spa- 
tial and time-reversal symmetries are exploited to restrict 
the Brillouin-zone (BZ) summations to irreducible sets of 
k-points. The nonlocal potential is calculated in the ba- 
sis of single-particle eigenstates, which allows to truncate 
the matrix at a certain number of bands. The divergence 
of the Coulomb potential in the BZ center is treated 
analytically instead of using dense k-point sets around 
k = 0. A nested density convergence scheme greatly re- 
duces the number of iterations in the self-consistent field 
cycle. Finally, by a suitable transformation of the MPB 
the Coulomb matrix becomes sparse, which speeds up 
the matrix-vector operations considerably. 

This transformation relies on the analytic properties 
of the bare Coulomb potential. Any other potential, in 
particular the screened Coulomb interaction, will not lead 
to a sparse matrix representation, though. Furthermore, 
in contrast to Gaussian or plane-wave basis sets, a direct 
evaluation of the screened Coulomb matrix, in the same 
way as for the bare Coulomb matrix,*^ is cumbersome 
in the MPB. Therefore, we incorporate the screening, 
after calculating the bare nonlocal exchange potential, 
in a separate step, which produces hardly any overhead. 
In this way, the simple analytic properties of the bare 
Coulomb potential as well as the sparsity of the Coulomb 
matrix are retained and can be taken advantage of. 

We validate our implementation by comparing results 
for prototypical semiconductors and insulators with re- 
sults from the literature. Then, we calculate the ground- 
state properties and the band structure of GdN. The 
band gap of GdN is controversially discussed in liter- 
ature. Results from LSDA-I-C/ calculations are incon- 
clusive. While the linearized muffin-tin orbital (LMTO) 
approach yields a narrow-gap semiconductor as ground 
state, "'^^ "^^ GdN exhibits a transition from a half- 
metallic to a semiconducting ground state under strain 
within the FLAPW method.C^l Two different solutions 
close in energy were obtained in an investigation using 
the hybrid functional BSLYP.fS^lBoth solutions were half- 
metallic, one in the majority spin channel, the other one 
in the minority spin channel. Our results show a transi- 
tion from a half-metallic ground state, which is similar to 
the energetically slightly less favorable solution of Ref. [231 
to a semiconductor under strain as in the LSDA-I-C/ cal- 
culation of Ref. TT, In their work, a large change of the 
lattice constant of more than 10% was necessary to ob- 
serve this transition. However, our calculations indicate 
that already small volume changes (~ 0.5%) are sufficient 
to observe this transition. 

The paper is organized as follows. In Sec. |llj we give a 
brief introduction to the theory of hybrid functionals. In 
Sec. |III[ we introduce the FLAPW method and describe 
our implementation of the HSE functional. In Sec. |IV[ 
we first compare results for prototypical semiconductors 
and insulators with values from the literature. Then, we 
present our findings for GdN in Sec. |Vj In Sec. |Vlj we 
draw our conclusions. 
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II. THEORY 

The construction of hybrid functionals as mixtures 
of local functionals with a nonlocal exchange term is 
motivated by the adiabatic-connection formula j^'^^'^^ ' in 
which the scaling of the electron-electron interaction es- 
tablishes a connection between the noninteracting Kohn- 
Sham system with the fully interacting one. In the 
weakly interacting limit, the formula becomes i dentic al 
to the HF exchange term which prompted Beck^^-'-^ to 
introduce a certain fraction a of HF exchange into the xc 
functional 



^HYB ^ 
XC xc 



a{E, 



HF 



(1) 



where E^^ denotes the local xc functional and E}^ its 
exchange part. E^^ is the HF exchange energy 



E. 



HF 



BZ occ. „ „ 
(T k,q n,7i' 

V^nUr)'P:'^{rH\v ~ r'|)^-q(r')^^k(r'), (2) 

evaluated with the Kohn-Sham wave functions V'5^k('^) 
spin a, wave vector k, and band index n. The sums 
over k and q run over the full Brillouin zone (BZ), n 
and n' sum over all occupied states, and v{r) — 1/r is 
the bare Coulomb potential. Here and in the following, 
atomic units are used unless stated otherwise. As the 
wave functions y'J^kl'") functionals of the effective po- 
tential, which in turn is a functional of the density, E^^ 
is a true functional of the density, too. 

Perdew et al^ deduced a mixing parameter a = 0.25 
by assuming a certain shape for the adiabatic-connection 
integrand. They proposed to use the Perdew-Burke- 
Ernzcrhof (PBE) functionaF* for the local part. The re- 
sulting functional 



E^ 



PBEO 
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HF 



E^ 
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(3) 



is nowadays referred to as PBEO. 

As the long-range part of the no nlocal HF term is 
cumbersome to calculate, Heyd et alW^ suggested to 
replace it again by a simple local functional. Later, 
it was demonstratecP^ that this leads to an improved 
description of the band gaps of semiconductors. Heyd 
et al.^^ used the error function erf(a;) and its comple- 
ment erfc(a::) = 1 — erf (x) to decompose the Coulomb in- 
teraction into a long-range (LR) and a short-range (SR) 
part 

r r 

where lo is an adjustable screening parameter. The HSE 
hybrid functional is thus given by 



PBE 



[E^^-^^ic.) - E^^'^^iu.)] , (5) 



where E^^'^^'{uj) corresponds to Eq. pi) with the bare 
Coulomb potential v{r) replaced by v^^r). E^^^'^^{uj) 



is the local functional for the SR exchange according to 
the decomposition given in Eq. (41). Its numerical treat- 
ment is discussed in Refs.[33]and[SlI Based on numerical 
fits to benchmark data sets of molecules, Heyd et al^ 
found an optimized value for the screening parameter 
uj — 0.15. In this work, we emp^lo^ the value of a; = 0.11, 
which was optimized for solids. 

Hybrid functionals are usually applied within the gen- 
eralized Kohn-Sham formalism,E21 which is based on a 
fictitious system of noninteracting electrons. These par- 
ticles move subject to a nonlocal potential that is defined 
in such a way that the electron density n(r) equals that 
of the real system. The nonlocal potential contains a lo- 
cal part that consists of the external potential created by 
the nuclear charges, the Hartree potential, i.e., the elec- 
trostatic potential produced by the total electron charge 
density, and a local contribution that derives from func- 
tional differentiation of the local parts of Eq. ^ and ^ 
for the PBEO and HSE functionals, respectively. The 
implementation of this local part of the xc potential re- 
quires only minor modifications of the DFT code, and we 
will focus on the nonlocal part in the following, which de- 
rives from the nonlocal exchange energy functional E^^ . 
Leaving out the scaling factor a, its matrix representation 
in the basis of Kohn-Sham eigenstates is given by 

BZ occ. 



q 7n 



For the HSE functional, the bare interaction v{r) would 
have to be replaced by the screened interaction v^^{r). 
Yet in practice, we first evaluate the nonlocal potential 
in the form of Eq. ^ and correct for the screening after- 
wards by subtracting v^^{r), as will be explained in the 
next section. 



III. IMPLEMENTATION 

A. Basis sets 

Our implementation of the HSE functional is based on 
the all-electron FLAPW method,'^^!^^!] which space is 
partitioned into non-overlapping atom-centered muffin- 
tin (MT) spheres and the remaining interstitial region 
(IR). The core states, which are confined to the spheres, 
are obtained from solving the fully relativistic Dirac 
equation. For the valence and conduction states we use 
a set of basis functions that are defined differently in 
the two regions of space: plane waves e'^*'+*^' '' with 
|k + G| < Gmax in the IR and linear combinations of 
ufp {r)Yim(j) in the MT spheres, where r is measured 
from the sphere center, uf^lr) are numerical functions 
defined on a radial grid, Y;m(r) are spherical harmonics 
with angular-momentum quantum numbers < I < Zmax 
and |m| < I, a labels the atom in the unit cell, and p 
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is an index for different radial functions. Gmax and Zmax 
are cutoff parameters. The linear combinations are such 
that the basis functions and their radial derivatives are 
continuous at the MT sphere boundaries. 

As in Ref. 1^ we evaluate the nonlocal potential, 
Eq. (|6|, with the help of an auxiliary basis {Mqj(r)} 
and its bi-orthogonal set {Mqj(r)}, where q is a Bloch 
vector and J is used to index these basis functions. By 
placing the completeness relation 



qj 



qJ 



(7) 



at both sides of v{r), the six-dimensional integral be- 
comes a sum over vector-matrix-vector products in the 
space of the MPB 

occ. BZ 

^x^fn'Ck) = - EEE(<kl<nk-qA^qj) 
m q / J 

xt;/j(q)(Mq^j^^k-ql¥'"'k>, (8) 

with the usual bra-ket notation {}\g) = J (i^r f* {r)g{r) 
and the Coulomb matrijp^ 

vij{ci) ^11 dVd3/M*,(r)z;(r,r')A/q,j(r'). (9) 



We need to evaluate the latter only once at the beginning 
of the self-consistency cycle. We note here again that the 
screening is accounted for in a separate step [cf. Eq. ( 11 )]. 

As Eq. ([s]) indicates, the auxiliary basis must be suffi- 
ciently complete in the space of wave-function products. 
We therefore construct this basis directly from products 
of LAPW basis functions. In the interstitial region, the 
auxiliary basis functions are plane waves with a new cut- 
off parameter G'^^^ and in the MT spheres, the basis con- 
sists of numerical functions of the form M1p{r) Y lm(^) 
with a cutoff imax- This so-called mixed product basis 
(MPB) can be systematically converged to represent the 
products of LAPW wave functions numerically exactly. 
In this way, the MPB is on the same level of accuracy as 
the all-ele ctron LAPW basis for the wave functions. It 
was 

showiJSIMl that G;„ax and L„ax can be well below 
their mathematically determined exact limit, i.e., 2G,nax 
and 2Zinax, and even below Gmax and /max, which makes 
the MPB a very efficient basi s. Furth er details about the 
MPB can be found elsewhere .EEHlll 



B. Implementation of the HSE functional 

It seems that the implementation of the HSE func- 
tional is now straightforward: the bare Coulomb poten- 
tial in Eq. ^ is replaced by the screened one and it is 
proceeded as in Ref. HI] for the case of the PBEO func- 
tional. However, in this way we would loose a very favor- 
able property of the bare Coulomb potential, namely, its 



multipole expansion, which makes it possible to render 
1 



1/r 
erfc(u;r)/r 
erf(a;r)/r 




FIG. 1: (Color online) Comparison of the bare (red solid line) 
and the screened (blue dashed line) Coulomb potentials with 
the difference between both (green dotted line). The differ- 
ence does not exhibit a divergence at r = and is a smooth 
function everywhere. Its Fourier transform converges very 
rapidly. 



the Coulomb matrix sparse by a simple unitary trans- 
formation of the MPB. The sparsity of w/j(q) speeds 
up the matrix-vector multiplications in Eq. ([s]) consid- 
erably. The screened Coulomb potential does not have 
this simple property. Furthermore, the direct evaluation 
of w/j(q) following the techniques of Ref. is not trans- 
ferable to the screened interaction. 

For these reasons, we evaluate Eq. ^ with the bare 
interaction as before and use a separate step to incorpo- 
rate the screening. This procedure is motivated by Fig.jl] 
which shows the bare and screened Coulomb potentials, 
v{r) and w^^(r), as well as the difference, v^^{r), as a 
function of the distance r , measured in units of Bohr radii 
flQ. While the two potentials diverge at r = 0, their dif- 
ference remains finite. It has a very smooth behavior for 
all distances and should thus be suitable to be described 
in reciprocal space. In fact, we find that only very few 
plane waves are needed to represent the difference accu- 
rately. 

To make use of the quickly converging Fourier expan- 
sion of the long-range potential v^^{r), we rewrite the 
total xc potential in the form 



T^HSE _ T/PBE 



PBE.SR 



NL 



NL.LRa 



(10) 



with the local xc potentials V^J^^ and ^pbe.sr ^-^^^^ 
rive from E^^^ and £'^^^'^^, respectively, and the po- 
tential terms in the parenthesis sum up to the SR nonlo- 
cal potential for HSE as 
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occ BZ 

K^n^'^'^^W = K^f^'(k) + ^^^{^nk\ ¥'mk-qXq+G) (Xq+G Xq+c) (Xq+GVmk-q| <'k) > (H) 

m q G 




^^mi''^^ 0^) ~ KTriri^C*)' ^''^ ^ functioii of the number of 
G vectors used for its construction. The convergence 
was studied for the case of bulk sihcon using a superceh 
with eight atoms. Machine precision is achieved with as 
few as 40 plane waves which would translate to ten for a 
primitive unit cell containing two atoms. This behavior 
is independent of the q point. 



10 20 30 

Number of G vectors 

FIG. 2: Root-mean-square (rms) deviation of the eigenvalues 



of the second term in Eq. (Ill from the fully converged result 
as a function of the number of plane waves in the Fourier 
transformation for the F point of silicon. We have used a 
supercell containing eight atoms. The 40 Fourier components 
would translate to ten in a calculation of the primitive unit 
cell with two atoms. 



where Xq+G(r) = e'^'^"''*^'' '" /VV is a plane wave normal- 
ized by the crystal volume V. We evaluate the Fourier 
transform of the wave-function products with the help of 
the MPB 



(<Pnkl<^mk-qXq+G> 



(<^«kl'^mk-q*V)(^VlXq+G) 



(12) 

where the first integrals on the right-hand side are calcu- 
lated routinely already for V^'^^^. The Fourier transform 
of v^^ is known analytically 



(Xq- 



-G \ v^^\ Xq+G') 



47r 



|q+G| 



Hq+G|V4w2 



(13) 

We note that any other form of the screened Coulomb in- 
teraction could easily be implemented at this stage. Since 
the matrix elements are diagonal in reciprocal space, the 
second term of Eq. (11) takes in practical terms negli- 



gible time to compute. From the fact that this function 
approaches zero very quickly with |q + G|, it is clear that 
the results are easily converged up to machine precision, 
even if the Fourier coefficient in Eq. (T2| falls off slowly 
with respect to |q+ G| because of the rapidly varying 
all-electron wave functions. Figure [2] shows the root- 
mean-square deviation of the eigenvalues of the matrix 



We note that the Fourier transform in Eq. ( 13 1 diverges 
as 1/ I q -I- G I in the limit q + G — >■ 0. The s ame d iver- 
gence is found for the bare Coulomb potential,'^^'^ such 
that the 1/ |q -I- G|^ terms cancel in the difference. The 
remainder is finite and is given by 



lim 

q+G^O 



47r 



q + G| 



1 



-|q+G|V4w' 



(14) 



We will later see that this nondivergent behavior of the 
screened interaction gives rise to a favorable k-point con- 
vergence. 

IV. VERIFICATION 

First, we present calculations for a prototypical set of 
semiconductors (C, Si, and GaAs) and insulators (MgO, 
NaCl, and Ar) and com pare the results with previous 
works from the literature! ^^ * ^^ * We focus in particular on 
direct and indirect band transitions. These are calcu- 
lated as the energy differences of the highest occupied 
and the lowest unoccupied eigenstates at the correspond- 
ing points in the BZ. We have taken the experimental 
lattice constants from Ref. [SS) In Fig. [3j we show the 
convergence of the band gap for silicon with the size of 
the k-point mesh. Within HSE this convergence is almost 
as fast as in PEE, whereas PBEO requires larger k-point 
meshes. This can be attributed to the nondivergent be- 
havior of the screened interaction at k = (s. Sec. Ill) 
and was already observed in Ref. [35l We find that an 
8x8x8 k-point mesh gives reliable HSE results for the 
band gap as well as for ground-state properties. 

In Table in we compare our results for the F — >■ F, F — >■ 
X, and F — > L transition energies with those obtained 
by the projector-augmented-wave (PAW) methocps' and 
experimental data. The band transitions are calculated 
for a set of materials at their experimental lattice con- 
stants with the functionals PBE and HSE. Overall, we 
observe excellent agreement between the calculated val- 
ues. In comparison to the experimental results, the HSE 
functional considerably improves on the PBE values. For 
semiconductors, the HSE transition energies are in very 
good agreement with experiment, while the larger gaps 
of insulators are still underestimated. 
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TABLE I: Kohn-Sham transition energies in eV obtained with the functionals PBE and HSE at experimental lattice constants 
compared with values from PAW calculations and experiment. An 8x8x8 k-point mesh was employed. 

This work PAW" Expt. 









PBE 


HSE 


PBE 


HSE 




GaAs 


r - 


^ r 


0.54 


1.43 


0.56 


1.45 


1.52,' 1.63" 




r - 


-> X 


1.47 


2.06 


1.46 


2.02 


1.90,' 2.01," 2.18" 




r - 


^ L 


1.01 


1.78 


1.02 


1.76 


1.74,' 1.84," 1.85" 


Si 


r - 


^ r 


2.56 


3.32 


2.57 


3.32 


3.05,"* 3.34-3.36," 3.4" 




r - 


-> X 


0.71 


1.29 


0.71 


1.29 


1.13," 1.25'' 




r - 


^ L 


1.54 


2.24 


1.54 


2.24 


2.06,^ 2.40" 


C 


r - 


^ r 


5.60 


6.98 


5.59 


6.97 


7.3' 




r - 


■> X 


4.75 


5.90 


4.76 


5.91 


— 




r - 


^ L 


8.46 


10.02 


8.46 


10.02 




MgO 


r - 


^ r 


4.77 


6.49 


4.75 


6.50 


7.7» 




r - 


-> X 


9.14 


10.86 


9.15 


10.92 






r - 


^ L 


7.93 


9.69 


7.91 


9.64 




NaCl 


r - 


^ r 


5.20 


6.57 


5.20 


6.55 


8.5'' 




r - 


-> X 


7.58 


9.05 


7.60 


8.95 






r - 


^ L 


7.30 


8.66 


7.32 


8.67 




Ar 


r - 


^ r 


8.70 


10.36 


8.68 


10.34 


14.2' 



"Reference [35l 
''Reference 1561 
'^Reference [57| 
''Reference j58| 
"Reference j59| 
■^Reference 
f Reference 1611 
^Reference 1621 
'Reference 1631 
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FIG. 3: (Color online) Convergence of the silicon band gap i5g 
for the functionals PBE (blue, dashed), PBEO (green, dotted), 
and HSE (red, solid) with respect to the k-point mesh (n x 
n X n). 



We compute the equilibrium lattice constants and bulk 
moduli by calculating total energies for different lat- 
tice constants and fitting the results to the Murnaghan 



TABLE H: Optimized lattice constants in A obtained with 
the PBE and the HSE functional. An 8x8x8 k-point mesh 
was employed. Results are compared to experimental results 
and calculations using the HSE functional within a paW^ 
and a GaussiarP^ method. 





This work 


PAW" 


Gaussian' 


Expt." 


Functional 


PBE 


HSE 


HSE 


HSE 




GaAs 


5.743 


5.660 


5.687 


5.710 


5.648 


Si 


5.472 


5.441 


5.435 


5.451 


5.430 


C 


3.571 


3.549 


3.549 


3.557 


3.567 


MgO 


4.265 


4.217 


4.210 


4.222 


4.207 


NaCl 


5.703 


5.627 


5.659 


5.645 


5.595 



"Reference 1351 
''Reference [52l 

"Experimental data taken from Ref. 1551 



we com- 



equation of state. In Table |n] and Table III 
pare the lattice constants and bulk moduli obtained with 
our implementation of the HSE functional with results 
from implementations based on plane-wave 
Gaussian basis sets.'^ The results of all three methods 
agree very well. For example the lattice constants calcu- 



7 



TABLE III: Bulk moduli in GPa obtained with the PBE and 
the HSE functional. An 8x8x8 k-point mesh was employed. 
Results are compared to experimental results and calculations 
using the HSE functional within a PAW method. 





This work 


PAW 


Expt.' 


Functional 


PBE 


HSE 


HSE 




GaAs 


64.5 


79.2 


70.9 


75.6 


Si 


88.9 


98.0 


97.7 


99.2 


c 


433 


467 


467 


443 


MgO 


153 


177 


169 


165 


NaCl 


21.3 


28.8 


24.5 


26.6 



" Reference 1351 

''Experimental data taken from Ref. 1551 



TABLE IV: Numerical parameters used for the calculation of 
GdN. 

parameter value 

k-point mesh 8x8x8 

mufSn-tin radii -RMT(Gd) = 2.33ao -Rmt(N) = 1.95ao 

plane- wave cutoffs Gmnx — ^.Qag^ G'^^^ — S.Gag^ 

angular-momentum /max(Gd) — 12 /max(N) — 10 

cutoffs Lmax(Gd) = 6 Lmax(N) = 4 

local orbitals" Gd: 5s,5p;7s,7p,6d,5f 

N: 3s,3p,4d,5f 
number of bands 200 

" Reference 1651 

lated in the FLAPW and PAW methods differ by less 
than 3 pm. Except for diamond, the HSE functional 
yields lattice constants and bulk moduli in much bet- 
ter agreement with experiment than the semilocal PBE 
functional, which tends to overestimate the former and 
underestimate the latter. 



V. GADOLINIUM NITRIDE 
A. Computational setup 

GdN crystallizes in the rocksalt structure, with a room- 
temperature lattice constant of acdN = 4. 
valence band of this material consists of the N 2s, 2p and 
the Gd 4/ states. The 4/ states are only half-occupied. 
The conduction band is formed by the Gd 5d and 6s 
states as well as the 4/ states in the minority channel. 

We determine the numerical cutoff parameters for the 
calculations in such a way that the difference between 
the total energies calculated at the experimental lattice 
constant, acdN, and at l.OlacdN changes by less than 
1 meV upon increasing the parameters. In Table [TV] we 
list the parameters used for the GdN unit cell (consisting 
of two atoms). In particular, we converged the k-point 



sampling, the size of the FLAPW basis, the size of the 
MPB, and the number of local orbitals. The latter are 
additional basis functions that are used to describe semi- 
core states'^ or to eliminate the linearization error .'^^21111 
In the following, we compare our theoretical HSE re- 
sults for the lattice constant, bulk modulus, band gaps, 
and magnetic moment with previous calculations and ex- 
periment. In order to compare our band structure results 
obtained at OK with the experimental results obtained 
at room-temperature, where GdN is in the paramagnetic 
state, we follow the idea that there is no difference be- 
tween the ferromagnetic and paramagnetic state for the 
exchange splitting and the large magnetic moments of 
the 4/ electrons. In the paramagnetic state we rather 
assume local magnetic / moments that fluctuate in di- 
rection with an overall zero magnetization. Thus, the 
magnetic polarization of the N states disappears. Disper- 
sive valence and conduction electrons that exhibit a large 
group velocity feel at any moment in time a small ran- 
dom potential landscape due to the exchange potentials 
of Gd / moments pointing in random directions. Follow- 
ing Ref. 46 this can be approximated assuming that in 
the paramagnetic phase each of these s, p, and d valence 
and conduction states characterized by a k-point band 
index are obtained by the averages of the correspond- 
ing spin-up and spin-down energies in the ferromagnetic 
phase. 

Furthermore, from the total-energy differences we de- 
rive the exchange coupling constants for a Heisenberg 
spin Hamiltonian and determine the Curie temperature. 
This gives us a measure for the quality of energy dif- 
ferences that can be expected from the HSE functional 
between different magnetic states. 



B. Structural and electronic properties 

We start our investigation of GdN by evaluating its 
structural and electronic properties and c omparing them 
to some of the available experimenta P" * ^^ -^^^ and theoret- 
ical data, obtained with the hybrid B3LYP functional^ 
and within the LSDA+C/ approach.^'' The compari- 
son is shown in Table |V] We note that our parameter-free 
HSE calculations yield a lattice parameter of 4.967 A in 
very close agreement to the experiment, while B3LYP 
overestimates the value by ~ 2% and the lattice con- 
stant in the LSDA-|-[/ method depends on the choice of 
the parameter U. Thermal expansion could account for 
the remaining difference to the experimental lattice pa- 
rameter that was determined at the room temperature 
(whereas the theoretical result corresponds to OK). The 
thermal expansion coefficient of GdN is unknown so far. 
Assuming linear expansion between K and room tem- 
perature (293 K) with the coefficient of isostructural and 
isovalent EuO (« 13 x 10~^ K"^),'^^ one would extrapo- 
late to 4.969 A at K, which is, indeed, very close to our 
optimized lattice constant. 

Next, we turn to the electronic structure. In Fig. [4] 
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TABLE V: Comparison of our HSE results for GdN with those from LSDA+C/ and B3LYP calculations and experiment. The 
theoretical results are given for the optimized lattice constant, unless stated otherwise. 







HSE LSDA+f/' LSDA+(7" 




T3 OT \/D c 

r>6\-i\V 


Expt. 


^ 

Lattice constant (A) 


(a no o \ 

(4.988) 


4.967 


4.92 


5.08 


(4.988) 


5.05 


4.988^ 


Bulk modulus (GPa) 




164 




150 




159 


192^ 


Magnetic moment (/ib) 


6.99 


6.99 




6.93" 




7.0 


6.88'' 


Direct gap at X (eV) (T < Tc) 


0.90 


0.85 


-0.16' 




0.91 


1.18' 


0.90^ 


Direct gap at X (eV) (T > Tc) 


1.17 


1.11 


0.10' 


0.98" 


1.30 


1.77' 


1.31^ 


Indirect gap L^X (eV) (T < Tc) 


0.01 


-0.06 


-0.45' 


0.14" 


0.43 


0.72' 




Indirect gap r-5>X (eV) (T > Tc) 


0.90 


0.85 


-0.13' 


0.69" 


0.98 


1.47' 




Position of majority 4/ peak (eV)* 


-6.00 


-6.00 


-7.8 


-8.1"'' 




-6.3' 


-7.8' 


Position of minority 4/ peak (eV)* 


6.05 


6.05 


6.6 


5.0"'' 


4.8' 


5.5' 


5.5 - 6.1*" 



"At the experimental lattice constant of 4.988 A. 

''Reference 1 14l U optimized for Gd bulk (Ref. I66I I. 

'^Reference 46 V chosen to reproduce the experimental direct gap 
of paramagnetic GdN (Ref. 67). 

''Reference 68 V chosen to reproduce the experimental direct gap 
of paramagnetic GdN. 

■^Insulating solution of Ref. 1231 

^At room temperature; Ref. 1691 

f Reference 70 

''Reference |10| 

'Extracted from the band structure. 
'Reference 1681 

*^Relative to the top of the valence band. 
'Reference 1111 

"Reference [7T| measured for GdX (X = P, As, Sb, and Bi). 




FIG. 4: (Color online) Band structure and density of states (DOS; in states per eV) of GdN at the experimental lattice constant. 
The majority and minority bands are plotted as solid and dotted lines, respectively. The orbital-resolved DOS is shown on the 
left for majority and on the right for minority states. The solid blue line shows the Gd 4/ states, the red dashed line the Gd 
5d states, and the green dotted line the N 2p states. 
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we show the spin-resolved band structure and the spin- 
and orbitally resolved density of states calculated at the 
experimental lattice constant (4.988 A). At about — 6eV 
and +6 eV we find the localized majority and minority 
Gd 4/ state, respectively. In the vicinity of the Fermi 
energy, GdN exhibits a truly interesting electronic struc- 
ture. Of particular interest are a direct and an indirect 
band gap, discussed extensively in the literature. The 
direct band gap accessible by optical measurements is 
located at the X point and amounts to 0.9 eV for the ma- 
jority spin channel and 1.5 eV for the minority spin chan- 
nel. The valence and the conduction states are separated 
by an indirect band gap (F— >-X): In the minority-spin 
states, this is a robust band gap of 1.5 eV, while in the 
majority-spin states this gap is tiny, only 0.01 eV. Thus, 
we observe that GdN is in a narrow-gap semiconducting 
ground state with an almost vanishing indirect band gap 
(F— >X) in the majority-spin direction, between the N 2p 
states in the valence and the Gd 5c? states in the con- 
duction band. This explains the different experimental 
reports disclosing GdN as a low carrier semimetal or an 
insulator depending on small changes of the experimental 
circumstances. 

Upon decreasing the lattice constant isotropically 
by just 0.02 A to the theoretically optimized value of 
4.967 A, we observe a transition to a half-metallic state 
or, more precisely, to a semimetallic state just for the ma- 
jority states: a small portion of the N 2p states at the F 
point remains unoccupied, while the Gd 5d band becomes 
partially occupied at the X point. If we define the "band 
gap" as the difference of these states, we formally get its 
value to be negative (see the band-transition energies in 
Table [v|). The described transition from a half- metallic to 
semiconducting state under isotropic strain was also ob- 
served by Duan et a/. however at a much larger lattice 
constant of 5.63 A. Our results suggest that, since the 
semiconductor to half-metal transition occurs so close to 
the equilibrium lattice parameter, the growth conditions, 
which influence the material properties such as the con- 
centration of possible N vacancies or strain in the system, 
or the lattice expansion upon temperature changes play 
a decisive role in the transport properties of GdN. 

Interesting physics can be expected also upon doping. 
If we n-dope into the conduction band or p-dope into 
the valence band, e.g. by use of Eu, at a concentration 
where electrons or holes populate only majority states, 
we can obtain significant charge currents with 100% spin 
polarization. In the paramagnetic phase, minority and 
majority states converge to spin-degenerate valence and 
conduction states, and any spin polarization of the charge 
current disappears. At the same time, the band gap 
opens in the paramagnetic state, partly because of the 
thermal expansion and partly because of the averaged- 
out exchange potential felt by the electrons. This, and a 
possible coupling of the conduction electrons or holes to 
the fluctuating 4/ momentJ^ may change the conduc- 
tivity by orders of magnitude when passing through the 
Gurie temperature (Tc). 



In Table [V] we list the band-transition energies calcu- 
lated with our method at the experimental and at the op- 
timized lattice constant for the ferromagnetic (T < Tq) 
and the paramagnetic {T > Tq) state. We find invariably 
larger band gaps for the paramagnetic state with no tran- 
sition to a metallic state nearby. Our calculated band- 
transition energies compare well with the experimental 
data where available. The small indirect majority spin 
gap between F and X of 0.01 eV compares well with esti- 
mates of 0.05 eV from Chantis et o/.'^ obtained using the 
quasiparticle self-consistent GW method (QSGW) com- 
bined with their empirical rule to estimate band gaps in 
semiconductors. The transition energies obtained within 
the LSDA-I-J7 method depend strongly on the choice of 
the parameter U. At the experimental lattice constant, 
we find similar transition ener gies as in the works of Lar- 
son et al^ and Trodahl et a/.p^ where in the latter the 
parameter U applied to the 4/ electrons was chosen to 
agree with the differences of the binding energies of the 
occupied and unoccupied 4/ level as measured with the 
X-ray photoemission and inverse photoemission, respec- 
tively, in the Gd pnictides and the U applied to Gd d 
states was chosen to reproduce the direct experimental 
band gap in the paramagnetic phase. In this way, the 
redshift of the direct band gap of 0.41 eV going from the 
paramagnetic state to the ferromagnetic state is perfectly 
reproduced, while our parameter-free calculation gives 
0.27 eV. The crossing of the conduction and the valence 
band at the X point was obtained with smaller values 
of whereas we find a direct band gap at X even 
below the Curie temperature Tq- The calculation with 
B3LYF^— yields three different solutions, where only the 
insulating one is similar to our result. The band gaps are 
significantly larger which may be attributed to the larger 
optimized lattice constant. 



The binding energy of the Gd 4/ majority band is also 
improved in the HSE scheme: the partial compensation 
of the self-interaction error leads to a pronounced shift of 
the localized 4/ states to larger binding energies. Calcu- 
lation with the PBE functional yields a much too shal- 
low / majority band, located at 3.1 eV below the Fermi 
energy; in HSE this band appears at a binding energy 
of 6.0 eV, much closer to its experimentally measured 
position at 7.8 eV.'^i' Furthermore, we note a very good 
agreement with the insulating B3LYP result, where the 
position of the 4/ peak is found at 6.3eV.'23 However, 
the agreement with experiment is not perfect. For the 
/ systems a stronger mixing of nonlocal exchange would 
probably give a better result. As compared to PBE re- 
sults, the unoccupied 4/ minority states shift towards 
higher energies (6.05 eV), which is cons istent w ith previ- 
ous LSDA+U and B3LYP calculations.'SIllllSI ^he posi- 
tion of the unoccupied 4/ states agrees well with typi- 
cal experimental results for the gadolinium pnictides be- 
tween 5.5 eV and 6.1 eV obtained with inverse photoe- 
mission spectroscopy.!^ 
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C. Magnetic order and critical temperature 

The ground state of unconstrained bulk GdN is ferro- 
magnetic (FM), with a Curie temperature of 58 Kl^and a 
magnetic moment of 6.88 /ie per Gd atonPSl determined 
from the saturation magnetization at 1.2 K. The total 
calculated magnetic moment is 7 /ib per formula unit of 
which 6.99 /iB comes from the Gd muffin-tin sphere. Our 
results are in good agreement with these observations 
(see the comparison of the theoretical and experimental 
values of the magnetic moment in Table |v|) . We confirm 
the magnetic order and determine the critical temper- 
ature by mapping the total energies obtained from our 
HSE calculations for several magnetic structures onto the 
classical Heisenberg Hamiltonian 



i y i=nn j=nnn J 



(15) 



including the nearest neighbors (nn), and the next near- 
est neighbors (nnn) interaction where Ji and J2 are the 
respective coupling constants with normalized spin vec- 
tors Si and Sj. Positive and negative values of J fa- 
vor parallel and antiparallel spin alignment, respectively. 
The coupling constants are extracted from the differences 
of the total energies of the FM, and two types of antiferro- 
magnetic (AFM) configurations characterized by planes 
of ferromagnetically ordered moments that are antiferro- 
magnetically stacked along the crystallographic [001] or 
[111] directions (AFM-I and AFM-II, respectively) .E^l For 
the calculation of the AFM-I (AFM-II) phase we use a 
tetragonal 1x1x2 (trigonal 0. x \/2 x 0) unit cell 
containing two Gd atoms and calculate the FM state in 
the same unit cells, in order to guarantee reliable total 
energy differences. All the calculations are performed at 
the experimental lattice constant. 
From the expressions 



AEu 



Eafm,i - 
Eafmji 



Efm,i = 8J1 
- EpMJI — 6J1 



6J2 



(16) 
(17) 



the Heisenberg coupling constants Ji and J2 are easily 
obtained. We list them in Table |VI[ along with their 
values calculated in previous studies using an LSDA-I-J7 
method within an FLAPW^'' and an LMTO basis set.^^ 
Both coupling constants are positive, confirming the fer- 
romagnetic nature of the ground state. 

We use two approaches to estimate the Gurie temper- 
ature. In the mean-field approximation (MFA) 



MFA 



1 

3h 



(12Ji + 6J2) , 



(18) 



we obtain a TJ^^^ = 55 K, very close to the experimental 
value of 58 K.^^^It is known, however, that the mean-field 
theory overestimates the Gurie temperature. For com- 
parison, we have also calculated the critical temperature 
by employing the random phase approximation (RPA) 



TABLE VI: Differences of total energies for different magnetic 
configurations [Eqs. (16 1 and (17 1], the Heisenberg coupling 
constants, and tfie corresponding Curie temperatures within 
the mean-field approximation [Eq. ( 18 1] and random phase 
approximation [Eq. (19|]. Energies and coupling constants 
are given in meV and the Curie temperatures in K. 





AEi 


AEu 


Ji 


J2 


rriMFA 


rriRPA 


This work 


8.8 


7.6 


1.09 


0.17 


55 


42 


Duan et al.'^ 


6.7 


4.2 


0.84 


-0.14 


36 


26 


Mitra et al} 


3.4 


0.4 


0.42 


-0.36 


11 


5 



"Reference [Til 
'Reference [12] 



as described in Refs. [751 and which is known to give 
results close to Monte Garlo solution: 



rriRPA 



1 

3A:b 



d^^ 



1 



Bz J(0) - J(q) 



(19) 



where we evaluate the integral on a discrete mesh of q 
points within the Brillouin zone. J(q) is the Fourier 
transform of the exchange coupling constants defined as 



J(q) = E Ji e'l-^"" 



(20) 



where Rnn and Rnnn are the positions of the nearest and 
the next nearest neighbors, respectively. The resulting 
T^PA = 42 K is roughly 30% smaller than the mean-field 
estimate. We consider these results as a sophisticated 
theoretical estimation of the Gurie temperature that goes 
along with a few uncertainties, some of which are difficult 
to assess, such as the quality of HSE being an approx- 
imation to the true but unknown exchange and corre- 
lation functional and to a much lesser extent the adia- 
batic approximation inherent in applying the Heisenberg 
model. Easier to assess are technically induced error es- 



timates: (i) AE in Eqs. (16 1 and (17) are converged to 



about 1 meV, which translates to an uncertainty of 3K. 
(ii) Base d on Monte Garlo calculations with two (as given 
in Table VI I and three nearest neighbors employing cou- 
pling constants published by Duan et aZ.,!^ which in both 
cases lead to the same Gurie temperature of 28 K, we es- 
timate that the neglect of exchange interactions beyond 
next nearest neighbors leads to a maximum uncertainty 
of 1 K. (iii) Employing our coupling constants (as given 
in Table fVll), we find that the RPA result of 42 K ap- 



proximates the numerically precise determination of the 
Gurie temperature within the Heisenberg model obtained 
by Monte Garlo, 45 K, by 3K. 

With these error estimates in mind, we compare our 
values to results of experimental studies, e.g. Tc — 68, 
69, 58 or 37 K as reported by Granville et aL^Khazen 
et al.^ Leuenberger et aZ.pl and Yoshitomi et al.^ 
which vary in value also depending on film thickness, 
strain, grain size, stoichiometry, and N vacancies. ^^ ' ^ 
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We conclude that our results are in very good agree- 
ment with the experimental situation. Comparing our 
results to other theoretical values exhibited in Table IVll 
we note that our coupling constants Ji and J2 obtained 
with HSE are significantly higher which gives rise to a 
higher Curie temperature, in agreement with experiment. 
We observe that the increase of the coupling constants 
goes along with an increase of the Gd 4/ moment in the 
muffin-tin sphere by 90m/ZB from 6.78 /xb to 6.87/LtB, a 
decrease of the Gd 5d moment by 20m/iB from 90m/XB 
to 70m/iB: and an increase of the N 2p moment, which 
is aligned antiparallel to the Gd 4/ moment, by 20m^B, 
from — lOOm/iB to — 120m/XB, when HSE is compared to 
the PBE functional. The precise understanding of the 
relationship between the change of the moments and the 
coupling constants requires additional analysis that goes 
beyond the scope of the paper. 

VI. CONCLUSION 

In this work, we have presented an implementation of 
the HSE hybrid functional, which contains a nonlocal 
screened exchange potential, within the FLAPW method 
as implemented in the Fleur code.'^ The calculation of 
the nonlocal exchange potential is realized by projecting 
the wave- function products onto the mixed product basis, 
reducing the six-dimensional integrations over the non- 
local interaction potential to vector-matrix-vector prod- 
ucts, where the matrix must be calculated only once at 
the beginning of the self-consistent-field cycle. 

We employ a sparse-matrix technique^ to evaluate 
the vector-matrix-vector products and incorporate the 
screening, i.e., the long-range part of the potential, in a 
separate step, where we exploit its fast converging Fourier 
series. This procedure allows constructing the nonlocal 
HSE potential from PBEO up to machine precision at 
a negligible computational cost. We note that this ap- 
proach is not restricted to the error function used in the 
HSE functional. In fact, our approach is quite general, it 
can be easily applied to arbitrarily screened interaction 
potentials. 

The results for lattice constants and band-transition 



energies obtained within our method show excellent 
agreement with previous results obtained with the 
FAWl^ and Gaussian-basecP^ methods. We have con- 
firmed the finding of Paier et al!^ that the k-point con- 
vergence within HSE is comparable to the conventional 
local FEE functional, whereas in FBEO much larger k- 
point meshes are necessary. 

In addition, we have calculated the properties of the 
rare-earth compound GdN. There is an ongoing discus- 
sion whether the ground state is insulating or metallic. 
In fact, within the HSE functional the ground state is 
very close to a phase transition: we observe a tiny in- 
direct band gap at the experimental lattice constant at 
room temperature, which vanishes at the theoretically 
optimized K lattice constant - the compound becomes 
half-metallic. The experimentally known band transi- 
tions are in good agreement with our theoretical results. 
Furthermore, we have calculated the coupling constants 
for the Heisenberg spin Hamiltonian from total-energy 
differences of ferromagnetic and antiferromagnetic con- 
figurations. The resulting Curie temperature of 42 K 
evaluated in the random-phase approximation is in good 
agreement with the experimental value of 58 K and gives 
confidence in the energetics obtained by HSE for different 
magnetic phases. From this we conclude that the HSE 
functional has the potential to describe the properties of 
rare-earth chalcogenides without the need for employing 
a Hubbard U parameter. We encourage the community 
to make use of the potential of the HSE functional to ex- 
plore the more subtle properties of the rare-earth chalco- 
genides such as the physics due to strain, dopands, or 
heterostructures. 
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